In [1]:
from IPython.core.display import HTML
css_file = '../ipython_notebook_styles/ngcmstyle.css'
HTML(open(css_file, "r").read())
Out[1]:
The aim of this computer laboratory is to gain a practical insight into the fundamental steps required to solve a partial differential equation (PDE) over a simple geometric domain using a finite element approach. In doing so, you will learn:
The symbolic‐numeric environment chosen for this laboratory is FEniCS which is a powerful and user‐friendly tool for solving PDEs. However, FEniCS is not compatible with Anaconda on a Mac it would seem; to use it, you must temporarily comment out the line setting up Anaconda in your .bash_profile. To get the notebook working together with FEniCS you need to
.bash_profilesudo easy_install pipsudo pip install ipython[all] --upgradeThis will be an old, but functional, version of the notebook that can be used when not using conda.
The Poisson equation is encountered in a wide range of domains in physics from heat conduction, electrostatics, diffusion and elasticity. The Poisson problem is defined as follows:
$$ \begin{align} - {\nabla ^2}u(x) &= f(x) & x & \text{ in } \Omega \\ u(x) &= {u_0}(x) & x &\text{ on } \partial \Omega \end{align} $$Here, $u(x)$ is the unknown function, $f$ and ${u_0}$ are prescribed functions, ${\nabla^2} = \Delta $ is the Laplace operator, $\Omega $ is the spatial domain and $\partial \Omega $ is its boundary.
A time-independent or stationary PDE like the Poisson equation combined with boundary conditions such as those above constitute a boundary-value problem.
Solving a physical problem with FEnics involves the following workflow:
Write a Python programme encoding the :
Add statements to the Python programme for:
The core of the recipe for turning a PDE into a variational problem is to multiply the PDE by a function $\upsilon $, integrate the resulting equation over $\Omega $, and perform integration by parts of terms with second-order derivatives.
The function $\upsilon$ which multiplies the PDE is in the mathematical finite element literature called a test function. The unknown function $u$ to be approximated is referred to as a trial function.
The terms test and trial function are used in FEniCS programs too.
The test function $\upsilon$ is required to vanish on the parts of the boundary where $u$ is known
Suitable function spaces must be specified for the test and trial functions. For standard PDEs arising in physics and mechanics such spaces are well known.
$$ \begin{equation} - \int_\Omega \nabla^2 u(x)\upsilon \, dx = \int_\Omega f(x)\upsilon \, dx \end{equation} $$Apply integration by parts to the PDE $ - {\nabla ^2}u(x) = f(x)$ with the divergence theorem (relation between the flow (flux) of a vector field though a surface and the behavior inside the surface) to recast it into its variational form. You will also take into account the conditions about the test function $\upsilon $ on the boundary where $u$ is prescribed.
You have now derived the variational formulation of the Poisson problem. This variational form is what is commonly called a weak form of the boundary-value problem ("weak" because of the reduced continuity requirements placed on $u$). It contains the basic equation $ - {\nabla ^2}u(x) = f(x)$ and the condition $u = {u_0} \text{ on }\partial \Omega $.
This continuous variational problem must be discretised to be solved using the finite element method. The continuous version of the unknown function $u$ will be denoted by ${u_e}$ (the index $e$ stands for "exact") to make the distinction with the discretised solution of the problem that is typically denoted by ${u_h}$.
NB: In the context of FEnics, $u$ will be meant to represent the discretised approximate solution of the weak form.
For a linear weak form like the one just established it proves convenient to introduce the following notation:
$$ \begin{equation} a(u,\upsilon ) = \int_{\Omega} {\nabla u \cdot \nabla \upsilon \, dx} \end{equation} $$$$ \begin{equation} L(\upsilon ) = \int_{\Omega} {f\upsilon \, dx} \end{equation} $$Where $a$ and $L$ are respectively a bilinear and linear operator.
The weak form of the Poisson problem (or any linear weak form) can therefore be expressed as:
$$ \begin{equation} a(u,\upsilon ) = L(\upsilon ) \end{equation} $$For every linear problem to be solved:
The explicit formulas for $a$ and $L$ are then coded in the programme.
Solving the discretised weak form consists in finding:
$$ \begin{equation} u \in V \text{ such that } a(u, \upsilon) = L(\upsilon) \quad \forall \upsilon \in V \end{equation} $$Note that
$$ \begin{equation} V = \{ v \in H^1(\Omega): v = u_0 \text{ on } \partial \Omega \} \end{equation} $$is the space of trial functions, and
$$ \begin{equation} \hat{V} = \{ v \in H^1(\Omega): v = 0 \text{ on } \partial \Omega \} \end{equation} $$is the space of test functions.
Specifying $V$ and $\hat{V}$ consists in choosing the mesh and the type of interpolation functions.
Here we are going to particularise the Poisson problem (reproduced here for convenience):
$$ \begin{align} - {\nabla ^2}u(x) &= f(x) & x & \text{ in } \Omega \\ u(x) &= {u_0}(x) & x &\text{ on } \partial \Omega \end{align} $$This decision means that $\{ {u_0},f,\Omega \} $ needs to be defined.
We choose a simple 2D domain: the unit square $\Omega = [0,1] \times [0,1]$ .
In order to compare the finite element solution to the exact solution we also a priori select a particular form for $u$. We choose:
$$ \begin{equation} {u_e}(x,y) = 1 + {x^2} + 2{y^2} \end{equation} $$If we inject equation the exact solution into Poisson problem we find that ${u_e}$ is solution if:
$$ \begin{equation} f(x,y) = - 6; \qquad {u_0}(x,y) = {u_e}(x,y) = 1 + {x^2} + 2{y^2} \end{equation} $$Using FEnics write an algorithm to solve the boundary-value problem described above.
In [ ]:
Here, we study the deflection $D(x,y)$ of an elastic circular membrane with radius $R$, subject to a localized perpendicular pressure force, modelled as a Gaussian function. The PDE representing this physical process is:
$$ \begin{equation} - T{\nabla ^2}D = p(x,y) \text{ in } \Omega = \{ (x,y)|\,{x^2} + {y^2} \le R\} \end{equation} $$where $p$ is the external pressure defined as:
$$ \begin{equation} p(x,y) = \frac{A}{{2\pi \sigma }}\exp \left[ { - \frac{1}{2}{{\left( {\frac{{x - {x_0}}}{\sigma }} \right)}^2} - \frac{1}{2}{{\left( {\frac{{y - {y_0}}}{\sigma }} \right)}^2}} \right] \end{equation} $$$T$ the constant tension in the membrane, $A$ the amplitude of the pressure, $({x_0},{y_0})$ the localisation of the Gaussian pressure function and $\sigma $ its width.
You will use the following values:
$$ \begin{align} T &= 10 \\ A &= 1 \\ R &= 0.3 \\ \theta &= 0.2 \\ {x_0} &= 0.6R\cos (\theta ) \\ {y_0} &= 0.6R\sin (\theta ) \\ \sigma &= 0.025 \end{align} $$
In [1]: